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Abstract 

We introduce a new variational method for the numerical homogenization of di- 
vergence form elliptic, parabolic and hyperbolic equations with arbitrary rough (L°°) 
coefficients. Our method does not rely on concepts of ergodicity or scale-separation 
but on compactness properties of the solution space and a new variational approach 
to homogenization. The approximation space is generated by an interpolation basis 
(over the nodes of a coarse mesh of resolution H) minimizing the norm of the 
source terms; its (pre-)computation involves minimizing 0{H~'^) quadratic (cell) 
problems on (super-)localized sub-domains of size 0{H\n{l/H)). The resulting lo- 
calized linear systems remain sparse and banded. The resulting interpolation basis 
functions are biharmonic (for c? < 3 and polyharmonic for d > 4) for the operator 
— div(aV-) and can be seen as a generalization of polyharmonic splines to differen- 
tial operators with arbitrary rough coefficients. The accuracy of the method {0{H) 
in energy norm) is established via the introduction of a new class of higher-order 
Poincare inequalities. The method bypasses (pre-)computations on the full domain 
and naturally generalizes to time dependent problems, it also provides a natural so- 
lution to the inverse problem of recovering the solution of a divergence form elliptic 
equation from a finite number of point measurements. 
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1 Introduction 



Consider the partial differential equation 

f-div (a{x)Vu{x)) = g{x) x e g e L^{n), a{x) = {uij e L°°(17)} 



u 



on 



(1.1) 



where is a bounded subset of with a piecewise smooth boundary (e.g., C^) and 
a is symmetric, uniformly elliptic on 0, and with entries in L°°(r2). It follows that the 
eigenvalues of a are uniformly bounded from below and above by two strictly positive 
constants, denoted by Amin(a) and Amax(a)- More precisely, for all ^ G M'^ and almost 
all X £ Q, 

Ami„(a)|el' < e'"a(x)e < A^ax(a)|ep. (1-2) 



In this paper, as in [55l [131 [56] , we are interested in the homogenization of (1.1) (and 
its parabolic and hyperbolic analogues in Section [s]), but not in the classical sense, 
(i.e., that of asymptotic analysis jl2j or that of G or i?- convergence ([52], |64| 



in which one considers a sequence of operators — div(aeV-) and seeks to characterize 



the limit of solutions. We are interested in the homogenization of (1.1) in the sense 



of "numerical homogenization," i.e., that of the approximation of the solution space of 



( 1.1 ) by a finite-dimensional space. This approximation is not based on concepts of scale 
separation and/or of ergodicity but on strong compactness properties, i.e., the fact that 
the unit ball of the solution space is strongly compactly embedded into 'Hq{0,) if source 
terms (g) are in (see [56] Appendix B]). 

In other words if ( 1.1 ) is "excited" by source terms g G then the solution space can 



be approximated by a finite dimensional sub-space (of ?^Q(r2)) and the question becomes 
"how to approximate" this solution space. The approximate solution space introduced 
in this paper is constructed as follows (see Section [2]): 

(1) Select a finite subset of points {xi)i^f>j- in Q, with := {1, . . . , N}. 

(2) Identify the nodal interpolation basis {(f)i)i(=j\/ as minimizers of | div(aV0) | subject 
to 4>i{xj) = 5ij (where is the Kronecker delta defined by 5i^i = 1 and (5jj = for 

(3) Identify the approximate solution space as the interpolation space generated by the 
elements (pi. 

Note that we construct here an interpolation basis with N elements, which inter- 
polates a given function using its values at the nodes {xi)i^j^ (analogously to e.g., the 
classical interpolation polynomials of degree A'^ — 1). 

Our motivation in identifying the interpolation elements through the minimization 
step (2) lies in observation that the origin of the compactness of the solution space lies 
in the higher integrability of 5. In Section [2] we will show that one of the properties of 



w 



the interpolation basis 4>i is that ^iWi(l)i minimizes | div(oViy)| over all functions 
interpolating the points Xi (i.e. such that w{xi) = wi). The error of the interpolation and 
the accuracy of the resulting finite element method is established in Section [3] through 
the introduction of a new class of higher order Poincare inequalities. The error analysis 



3 



is performed by identifying the points Xi as the nodes of a mesh Th of resolution H 
(although our method is meshless in its formulation such a mesh is needed to characterize 
the distribution of the points Xj, the distribution of those interpolation points can be 
increased near points of the domain where high accuracy is required). In the fully discrete 
formulation of our method (Section [o]), a will be chosen to be piecewise constant on a 
sub-mesh of Th of resolution h <^ H {in numerical applications where a is represented 
via constant values on a fine mesh of resolution h, H will be the resolution of a coarse 
mesh having interior nodes Xj). 

In Section [4] we show that our method also provides a natural solution to the (in- 
verse) problem of recovering u from the (partial) measurements of its (nodal) values 
at the sites Xj. Although we restrict, in this paper, our attention to d < 3, we 
show, in Section [5j how our method can be generalized to d > 4 by requiring that 
g £ L'^^{^) (for {d — l)/2 < m < d/2) and obtaining the elements 4>i as minimizers 
of I div(aV(/>) 1^™. The resulting elements are 2m-harmonic in the sense that they 

satisfy ( div(aV-))^'"0i = in il/{xj : j G Af}. In that sense our method can be seen 
as a polyharmonic formulation of homogenization and a generalization of polyharmonic 
splines |40^ [20] to PDEs with rough coefficients (see Section [6| the basis elements (f>i can 
be seen as "rough polyharmonic splines" ) . In Section [7] we show how the computation 
of the basis elements 0j can be localized to sub-domains i7j C and obtain a posteriori 
error estimates. Those a posteriori error estimates show that if the localized elements 
(/)^°^ decay at an exponential rate away from Xj (which we observe in our numerical ex- 
periments in Section [9]) then the accuracy of the method remains 0{H) in "^^^-norm for 
sub-domains Qi of size 0{Hln{l/H)^, which we refer to as the super-localization prop- 
erty. We defer the proof and statement of a priori error estimates to a sequel work. We 
give a short review of numerical homogenization in Section 10 We discuss and compare 
the accuracy and computational cost of our method in Subsection 10.3 The linear sys- 
tems to be solved for the identification of the (super-)localized basis ip^"^ are sparse and 
banded (by virtue of this property, the computational cost of our method is minimal) . 



2 Variational formulation and properties of the interpola- 
tion basis. 

2.1 Identification of the interpolation basis. 

Let V be the set of functions u G ^o(^^) such that div(aVu) G L'^{Q). Observe that V 



is the solution space of (1.1) for (7 G L^(0). By [65l|33], solutions of (1.1) are Holder 
continuous functions over Q, provided that their source terms g belong to L'^^'^'^'^ (Cl) for 
some £ > 0. It follows that if the dimension of the physical space is less than or equal to 
three (d < 3) then elements of V are Holder continuous functions over Q. For the sake 
of simplicity, we will restrict our presentation to d < 3 and show in Section [5] how the 
method and results of this paper can be generalized to d > 4. 
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Let ||.||y be the norm on V defined by 

\\u\\v := II div(aVn)||i2(Q). (2.1) 

It is easy to check that {V, \\.\\v) is a complete linear space; in particular, it is closed 
under the norm ||.||y. 

Let H > he a parameter determined by available computing power and desired 
accuracy. For simplicity, we assume that is a polygonal domain. Let Th be a coarse 
shape-regular tessellation of il. of resolution H, i.e. a tessellation (triangulation for d = 2) 
of 17 with tets (tetrahedra for d = 3, triangles for d = 2) of uniformly bounded aspect 
ratios, such that sup^-g-j-^ diam(T) < CH and such that the boundary of Th coincides 
with that of O. Write {xi, i = 1, . . . ,N) the non-boundary nodes (vertices) of Th and 
M the corresponding set with 

N := \M\ (2.2) 

For each i £ {1, . . . , N} define 

Vi:={(j)eV\ (pixi) = 1 and (l){xj) = 0, for j e {!,..., N} such that j / i} (2.3) 
and consider the the following optimization problem over T^: 



{Minimize | div(aV(/))| 
Subject to (f) ^ Vi. 



(2.4) 



We will now show that (2.4) is a well posed (strictly convex) quadratic optimization 



problem and we will identify its unique minimizer (pi. 

Remark 2.1. Observe that (pi depends only on a, and the locations of the points 
{xi)i(^j\f (the mesh Th is required only for error analysis). Note that we do not require 
the mesh Th to be quasi-uniform, i.e. we do not require its triangles to be quasi-uniform 
in size (the idea is that the density of points Xi may be increased in locations where higher 
accuracy is desired). 



Let G{x,y) be the Green's function of (1.1). Recall that G is symmetric and that it 
satisfies (for y £ 



\- div [a{x)VG{x,y)j = 6{x - y) x e 
\G{x,y) = for x G dQ, 

where 5{x — y) is the Dirac distribution centered at y. For i, j G {1, . . . , A^} define 

0jj := / G{x,Xi)G{x,Xj) dx. (2-6) 
Jn 
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Lemma 2.2. The integrals (2.6) are finite, well defined and for all i,j G {1, . . . ,N}, 
Qij is bounded by a constant depending only on Amin(a), Amax(a) and Q. The N x N 
matrix @ is symmetric positive definite. Furthermore for all I G M.^ , 



l^Ql = Ml2^n) (2-7) 



where v is the solution of 



— div ya{x)'Vv{x)j = ^f=i ljS{x — xj) for x £ 
t'(x) = for X G 



(2.8) 



Proof. By Theorem 1.1. of [39], for d> 3, the Green's function is bounded by 

G{x,y)<Kd\x-y\^-^ (2.9) 

where depends only on d, Amin(a) and Amax(o)- For d = 2, we have (see for instance 
|67| and references therein) 

Gix,y)<K,{l + ln^^^^) (2.10) 

\x — y\ 

where K2 depends only on Xmm{a), Amax(o) and ^l. For d = 1, the Green function is 

trivially bounded by a constant depending only on Ainin(a)) Amax(o) and 0. It follows 
from these observations that for d = 1,2,3 

/ {G{xi,y)fdy<K (2.11) 
Jn 

where the constant K is finite and depends only on Aniin(o), Amax(a) and fi. We deduce 
that the integrals in (2.6) (and hence Q) are well defined. 
Let / e M^, write 

N 

v{x):=^G{x,x^)li (2.12) 

i=l 



Observe that v is the solution of (2.8) and that ||^||^2(f^) = ^"^0^- Then, it follows that 
O is symmetric positive definite. Indeed if Q is not positive definite, then there would 
exist a non zero vector I G such that @l = 0. This would imply ||f ||l2(-q) = which 

is a contradiction since the equation — div ^a(x)Vw(x)^ = Ylf=i Ij^i^ — xj) has a non 
zero solution. □ 

Let P be the N x N symmetric definite positive matrix defined as the inverse of Q, 

i.e. 

P:=e-^ (2.13) 



Note that Lemma |2.2| (i.e. the well-posedness and invertibility of 0) guarantees the 
existence of P. 
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Define 



■r{x,y) := / G{x, z)G{z,y) dz. 



(2.14) 



Note that r is symmetric and (see proof of Lemma [2^ that for each y G Q, x — )• t{x, y) 
is a well defined element of V (in particular, it is Holder continuous on $7). Note also 
that for i^j G M we have T{xi,Xj) = 
(div(aV-)) in the sense that for y G 1^ 



Qij and that r is the fundamental solution of 



div ^aV( div(aVr(x, y))) j = 6{x — y) for x G 17 
T{x,y) = div (aVr(x,y)) =0 for x G d^. 



(2.15) 



Proposition 2.3. Vi is a non-empty closed afftne subspace of V . Problem (2.4) is a 
strictly convex optimization problem over Vi. The unique minimizer of (2.4) is 

N 



cpiix) := '^Pij t(x,; 



(2.16) 



Remark 2.4. It is important to note that in practical (numerical) applications each 
element (pi would be obtained by solving the quadratic optimization problem (2.4) rather 
than through the representation formula (2.16). 

Remark 2.5. Note that the tessellation Th is not required to construct the elements (j)i 
(it is only needed for error analysis). In that sense our method is meshless and depends 
only on the location of the nodes Xi ( and in practical applications the density of points 
can be adapted to the structure of a). More precisely, in the constructions below, we do 
not use the triangulation or any other tessellation, we just use the set of nodes. 

Proof. Let us first prove that G Vi. First observe that for all j G {1, . . . , N}, 

- div (^a(x)Vr(x, Xj)) = G(x, Xj) (2.17) 



and r(x, Xj 
Lemma 
that 



on 90. Noting that || div(a(x)Vr(x, Xj)) ||j^2(q) = ©jj we deduce from 



2.2 



that r(x,Xj) G V. We conclude from ( 2.16[ ) that 0j G V. Now observing 



G 



T(^Xi, Xj) 



we deduce (using the fact that P • G is the N x N identity matrix) that 

(piixj) = {P ■ Q)ij = 5i. 



(2.18) 



(2.19) 



We conclude that 4>i €Vi which implies that Vi is non empty (it is easy to check that it 
is a closed affine sub-space of V) . 

Now let us prove that problem (2.4) is a strictly convex optimization problem over 
Vi. Let v,w eVi such that v w. Write for A G [0, 1], 



/(A) := / dw{aV{v + X{w-v))) 



(2.20) 
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and we need to show that /(A) is a strictly convex function. Observing that 



/(A)= [ |div(aVw)f + 2A [ (div(aVv))(div(aV(u; - w))) + A^ [ \diviaV{v - w))f 
Jn Jn Jn 

(2.21) 

and noting that | div(aV(f — w))\ > (otherwise one would have v = w given that 
V and w are both equal to zero on cJfi) we deduce that / is strictly convex in A. We 



conclude that (see, for example, [28, pp. 35, Proposition 1.2]) that Problem (2.4) is a 
strictly convex optimization problem over Vi and that it admits a unique minimizer in 



Let us now prove that is the minimizer of (2.4). Let v Vi with v ^ Write 
I{v) := I div(aVf)p. We have 



I{v) = I{(t)i) + I{V - 00 + 2 J((^i, V - (t>^) 



(2.22) 



with 



J{<Pu v-^i)= I (div(aV0i))(div(aV(t; - <i)i))). (2.23) 
In 



Note that (as before) I{v — (pi) > ii v ^ (pi and using (2.17) and ( |2.16 ) we obtain that 

N 



div(a(x)V</.i(x)) := ^P,jG(a;,; 



which leads to 



N 



(2.24) 



(2.25) 



J{(j)i,v- (j)i) = '^Pij {v{xj) - Mxj)) = 

where in the last equality we have used v{xj) — <pi{xj) = for all j G {1, . . . , N}. It 
follows that I{v) > I{4'i) which implies that (j)i is the minimizer of (2.4). □ 



2.2 Variational properties of the interpolation basis 

Write Vq the subset of V defined by 

Vo:={vGV: vixi) = 0, Vi G {1, . . . , iV}} 
For u,v € V, define the (scalar) product {■,-) by 



{u,v) := / (div(aVn))(div(aV?;)) 
Jn 



(2.26) 



(2.27) 



The following theorem shows that the functions (f>i generate a linear space interpo- 



lating elements of V with minimal (^., .)-norm (see (2.29)). 
Theorem 2.6. It holds true that 
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The basis (f)i is orthorgonal to Vq with respect to the product (•,•), i-( 



0, \fie{l,...,N} and\fv G Vq 
• 'Wi(j)i is the unique minimizer of 

v2 



n 



(div(aVu;))^ 



(2.28) 



(2.29) 



over all w £ V such that w{xi) = Wi. 
For all i E {1, . . . , A^} and for all v £ V, 



N 



For all i,j S {1, . . . , A^}, 



(2.30) 



(2.31) 



Proof. Equation (2.28) trivially follows from (2.30). (2.30) is a direct consequence of 
(2.24). (2.30) leads to {4>i,v) = for w G Vq. To prove that w = Yli'^i^i is the unique 
minimizer of (2.29) over W := {w G V \ w{xi) = Wi for all i G {1, . . . ,-/V}} we simply 
observe that if w' is another element of W then w' — w belongs to Vq and for w' ^ w we 
have (from (2.28), using < w,w' — w >= 0) that 

(^w' , w''j = (^w, w'j + (^w' — w,w' — w'j > (^w, w'j . (2.32) 

ve have 

(2.33) 

□ 



For (2.31), using (2.30) and writing 1^ the N x N identity matrix we have 

N 



k=l 



Remark 2.7. It is easy to check that the orthogonality (2.28) implies that (pi solves 

div (^aV(div(aV(/)i))) =0 on n/{xj \ j £ M} 
(pi = div(aV(/)i) = on on (2-34) 

^(t>i{xj) = 6ij. 



However, (2.34) alone cannot be used to uniquely identify (pi as (2.16). This is due to 
the fact that the solution of (2.34) may not be unique. As a counter-example consider 
Q. = (—1, 1), xi = 0, A/" = 1 and a = I^. For a G M define 

jipaix) := '^x^ + |(q - l)x2 + ax + l on [-1,0] 
I ^P^(x) := ^x^ - |(a + l)x'^ + ax + I on [0, 1] 



(2.35) 
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then for any a £ M, is a solution (2.34). Using the variational formulation of (pi it 



is possible to show that to enfo rce t he uniqueness of the solution of (2.34) we also need 
to require that (see Proposition 7.6), for some c G M^, 



N 

div ^aV( div(aV(/)i))^ = ^^Cj5(x — xj) on Q 



(2.36) 



which in dimension one is equivalent to the continuity of div{a'V(j)i) across the coarse 
nodes {xj)j^j\f. Note that (2.36) implies that for all v £Vq 



j div(aV0i) div(aV?;) = j u div (^aV( div(aV</>i))) = 0. (2.37) 

3 From a Higher Order Poincare Inequality to the accu- 
racy of the interpolation basis 

In this section we introduce a new higher order Poincare inequahty and derive the 



accuracy of the interpolation space computed in (2.4). This new class can be thought 



of as a generalization of Sobolev inequalities for functions with scattered zeros (see [53] ) 
to operators with rough coefficients. 

3.1 A Higher Order Poincare Inequality 

We will first present the new Poincare inequality: 
Theorem 3.1. Let f G Vq. It holds true that 

l|V/k2(f,)<OT||div(aV/)||^,(^) (3.1) 



where the constant C depends only Amin(a) and on the (uniform) hound on the aspect 
ratios of the elements (tetrahedra for d = 3, triangles for d = 2) of Th- 
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Proof. Write \Vf\l := {VfY'aVf, 



|V/|^ = - / /(divaV/) (3.2) 



(3.3) 



= - E / /(divaV/) 

< E(/ /')'^'(/ |divaV/|Y/2 (3.4) 

<CY.H{I \Vf\'fl\j |divaV/|2)V2 (3.5) 

<CH{Y, l\Vf\')'/\Yl /jdivaV/lY/^ (3.6) 



<CF(/ \Vf\^f'\f |divaV/|2)V2 (3.7) 
where (3.1) follows. In the above, we have used the following inequaltiy, 

f f< CH^ [ |V/|2 (3.8) 
Jt Jt 

which is a standard property of piecewise linear interpolation |3H pp. 12, Proposition 
1.12]. Namely, \f - Xi/p < CH^ \Vf - VXi/p where Iif is the piecewise linear 
interpolation of /, and in this case Xi/ = 0. □ 

3.2 Accuracy of the interpolation basis 



Now let us use (3.1) to estimate the interpolation error. 



Corollary 3.2. Let u G V be the solution of (1.1) and ii™(x) := J2f=i'^{^i)4'i(,x)- H 
holds true that 

h-^'''\\nlin)<CH\\g\\^.^^^ (3.9) 

where the constant C depends only on Amin(a) and on the (uniform) bound on the aspect 
ratios of the elements (tetrahedra for d = 3, triangles for d = 2) ofTn- 



Proof. Observing that u — G Vo(f^) we deduce from Theorem 3.1 that 



/ \Viu-u''')\'^ <CH^ [ |5- Vn(j;i)divaV(/)i|2(ix (3.10) 
Jn Jn ^ 

<CH^i[ g^dx+ [ (Vn(xi)divaV(/>i)2(irr) (3.11) 
Jn Jn '~ 

<CH'^{ \ g'^dx+ I {d\YaVufdx) (3.12) 
Jn Jn 

< CH^gWl^^^^ (3.13) 
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Note that in the third inequahty, we have used the Lemma 2.6, which states that the 
hnear combination of (pi minimizes || • ||y among all the functions in V{0,) sharing the 
same values at the nodes {xj}. Note also that through this paper, for the sake of clarity, 
we we only keep track of the dependence of C and not its precise numerical value (we 
will, for instance, write 2C/Ainin(o) as C to avoid cumbersome expressions). □ 

3.3 Accuracy of the FEM with elements 0j 

The following theorem shows that the finite element method with elements (pi achieves 



the optimal (see |13| ) convergence rate 0{H) in T-l^ norm in its approximation of (1.1). 
Theorem 3.3. Let u be the solution of equation (1.1), and be the finite element so- 



lution of ( 1.1 ) over the linear space spanned by the elements {(pi, i = 1, . . . , A^} identified 



in (2.4). It holds true that 



\u — u 



Hi 



I "Hi 



lin)<CH\\g\\j^.^^^ (3.14) 



where the constant C depends only on \rmn{o), Amax(a) (^nd the (uniform) bound on the 
aspect ratios of the elements (tetrahedra for d = 3, triangles for d = 2) ofTn- 



Proof. The theorem is a straightforward consequence of Corollary |3.2| and of the fact 
that minimizes the (squared) distance (Vn — Vu^)^a{'S/u — Vu^) over the linear 
space spanned by the elements (pi. □ 



Remark 3.4. Recall that the convergence rate of a FEM applied to (1.1) with piecewise 



linear elements can be arbitrarily bad 11 Oj . Recall also that the convergence rate of 



Theorem 3.3 is optimal f73] / (this is related to the Kolmogorov n-widths [UB/ . 
which measure how accurately a given set of functions can be approximated by linear 
spaces of dimension n in a given norm). 

4 Recovering u from partial measurements 

In this section we will show that our method provides a natural solution to the inverse 



problem of recovering u from partial measurements. Consider the problem ( 1.1 ). In this 
inverse problem one is given the following (incomplete) information: a(x) is known, g is 
unknown but we know that Hs'llia^Q) < M, u is unknown but its values are known on a 
finite set of points {xj : j £ M} (through site measurements). The inverse problem is 
then to recover u (accurately in the Ti^ norm). A solution of this inverse problem can 



3.2 



be obtained by first computing the minimizers of (2.4) (or thei r loc alized version (7.3)) 
and approximate u with u™(x) := YliLi''^i^i)4'i{x)- Corollary 
bound the accuracy of the recovery by 



can then be used to 



\\u-u'^^.^^^^<CHM (4.1) 

where the constant C depends only on Amin(a) and on the (uniform) bound on the aspect 
ratios of the elements of a tessellation of the nodes Xi. Note that the method is meshless 
and that the tessellation is only required to estimate the constant C. 
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5 Generalization to <i > 4 

Although for the sake of simphcity we have restricted our presentation to d < 3, the 
method and results of this paper can be generalized to d > 4 by introducing the space 
V"^ defined as the set of functions u G 'Ho(r2) such that ( div(aV-))'"M G L^(0) where m 
is an integer m > 1 and (div(aV-))" u is the m-iterate of the operator (div(aV-)), i.e. 

(div(aV-))^n := (div(aV'u)) and ( div(aV-))'"u := ( div(aV-))™"^ ( div(aVu)) . (5.1) 
Introducing as the norm on defined by 



'dw(aV-)y'u . (5.2) 



The results of this paper can be generalized to d > 4 by choosing m such that (d — 1)/2 < 
m < ci/2, assuming that g € l?'"^(yt) and using the interpolation basis defined as the 
minimizer of 



(5.3) 



{Minimize HfPHv^ 
Subject to G VI 



where 

■= {</) G y"^ I ^{xi) = 1 and (\>{xj) = 0, for j G {1, . . . , A^} such that 3 / i}. (5.4) 



Note that the solutions of (5.3) are 2m-harmonic in the sense that they satisfy 



div(aV-))^™'</>(x) = for X 7^ Xy (5-5) 



Note that the solutions of (2.4) are bi-harmonic on rj/jxj : j G AA} 



For d = \ one can also obtain an accurate interpolation basis by minimizing Vcfya^cf) 



subject to the pointwise interpolation constraints used in (5.4). By doing so one would 
obtain basis elements (pi that are harmonic in Vt/{xj : j G A^} and recover a numerical 
homogenization method based on "harmonic coordinates" [71 155j. 

6 Rough Polyharmonic Splines 

6.1 The elements 0^ as rough polyharmonic splines 

The interpolation basis (pi constructed in this work can be seen as a generalization of 
polyharmonic splines to differential operators with rough coefficients. More precisely, as 
polyharmonic splines lead to an accurate interpolation of smooth functions (solutions of 
the laplace operator), the "rough polyharmonic splines" introduced in this paper lead 
to an accurate interpolation of solutions of differential operators with rough coefficients. 
Note also that as polyharmonic splines are m-harmonic in the sense that they satisfy 



A'^cp = away from the coarse nodes, the solutions of (5.3) are 2m-harmonic in the 
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sense that they satisfy ( div(aV-)) '"(/'(x) = away from those coarse nodes. This gener- 
alization is chahenging in several major ways due to the lack of regularity of coefficients. 
In particular, the Fourier analysis toolkit is no longer available and the simple task of 
enforcing clamped boundary conditions (trivial with smooth coefficients) becomes a sig- 
nificant difficulty with rough coefficients: consider, for instance, the problem of finding a 
function G Ti^ {^{0, 1)) that satisfies zero Neumann and Dirichlet boundary conditions 
on dB{0, 1) and such that div(aVi?^) G L^(i?(0, 1)). We will now give a short review of 
polyharmonic splines. 



6.2 Polyharmonic splines: a short review. 

Let X be a discrete (possibly finite) set of points of M"'. Let n be a real valued function 
defined on X. Polyharmonic splines have emerged through the problem of interpolating 
u (from its known values on X) to a real function (j) defined on M'^ (with (t){xj) = u{xj) for 
Xj € X). Let m be an integer strictly greater than d/2. The idea behind polyharmonic 
splines is to seek the most "regular" interpolation by selecting i?i> to be a minimizer of 
the semi-norm 

E -^'^Mfdxf (6.1) 

over all functions ip iii 'H"*(M'^) interpolating u on X^ where the parameters Cq- are 
positive coefficients usually |2Tl |62] chosen to be equal to ^ to ensure the rotational 
invariance of the semi-norm (those coefficients are also sometimes chosen to be equal 
to one [601 E])- Writing A the Laplace operator on M"^ (Acf) = Yfi=i 0) and A*" the 
m-iterate of </> (A^ = A and /^^(p = A(A'=-V)) the solutions of (lof satisfy A^c/) = 



on M. /X and are therefore called m-harmonic splines. A Polyharmonic spline of order 
m is therefore also commonly defined [13] as a tempered distribution (f) on that is 
2m — d — 1 continuously differentiable and such that /SJ^cj) = on / X. 

Polyharmonic splines can be represented via weighted sums of shifted fundamental 
solutions of A, (i.e. they can be written [20j (p = Ylx &x CjT{x — Xj) +pm-iix) where r 
is the fundamental solution of A"* (A™r(x) = S{x)) and Pm-i is a polynomial of degree 
at most m — 1). 

When X forms a regular lattice of M'^ the resulting splines are referred to as polyhar- 
monic cardinal splines. In this case (which has been extensively studied [631 145 ^ H6 l I47j) 
it can be shown |591 [60l |6T] that the basis elements allowing for the interpolation of u 
can be expressed as A^^r (where r is the fundamental solutions of A™ and A^ is the 
finite difference discretization of A on Z'^). 

Polyharmonic splines can be traced back to the seminal work of Harder and Des- 
marais jlQ] on the interpolation of functions of two variables based on the minimization 
of a quadratic functional corresponding to the bending energy of a thin plate (we refer 
to [T7] for a review) and the extension of this idea to higher dimensions in the seminal 
work of Duchon [201 l2Tl |22] ( building on earlier work by Atteia [6] ) . We also refer to 
the related work of Schoenberg {d = 1) on cardinal spline interpolation |63j . 
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7 Localization of the interpolation basis 
7.1 Introduction of the localized basis 

For each i G {1, . . . , N} let be an open subset of 0, containing Xi. Let AA(rij) be the 
set of indices corresponding to the interior nodes of the coarse mesh contained in $7j (i.e. 
M{Qi) ■= {j e {l,...,N}\xj e Qi}). Define 

V{n^) := {(P € nlin) \^ = OondniU (n/Qi) and dw{aV(j)) G L'^{n^)} (7.1) 

and 

Vii^i) :={(/) G y(Oi) I (/>(xi) = 1 and (t){xj) = 0, for j G 7V(Sli) such that j / (7.2) 
Consider the the following optimization problem over Vi{ni): 



{Minimize | div(aV(^)| 
Subject to G Vi{^li). 



(7.3) 



The proof of the following proposition is identical to that of Proposition 2.3 {Vi is 
simply replaced by Vi{Qi)). 

Proposition 7.1. Vi{il.i) is a non-empty closed convex affine subspace ofV{i}i) under 
the norm := Jq. (div(aVf)) . Problem (7.3) is a strictly convex (quadratic) 

optimization problem over Vi{Qi) with a unique minimizer denoted by 



Remark 7.2. As in Proposition \2.S\ the unique minimizer of (7.3) can be represented 
using the Green's function of the operator — div(aV-) with Dirichlet boundary condition 
on dili, i.e., we have for x £ Cli 



with 



•'°^(x,y):=/ G^''°^(x,z)G*>^(z,y)dy (7.5) 

JQ.i 

where G*'^°^(x,2/) is the Green's function of — dw{a'S/) in Qi with zero Dirichlet boundary 
condition, i.e. the solution of (for y E Qi) 

i-div(^a{x)VG^'^°%x,y)^=6{x-y) for x £ Qi; 
\G'>=(x,y) = for x G 9^, 

and, writing Ni := |AA(J7j)| the number of coarse nodes xj that are contained in 
pj,ioc ]\^. ^ positive definite symmetric matrix, defined over J\f{Qi), by = 

(^0*,ioc-j-i yjfiQj-g^ Q«,ioc X Ni positive definite symmetric matrix, defined over 

MiQi), by 

0l'!r-=/ G'''°'{x,Xk)G''''''{x,Xj)dx. (7.7) 
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Remark 7.3. It is important to note that in practical (numerical) applications each 
(localized) element cj}"^ is obtained by solving one (local, i.e. over Qi) quadratic op- 



timization problem (7.3) rather than through the representation formula (7.4) (which 
would require solving Ni elliptic problems in Vti corresponding to the values of the local 
Green's function). 

7.2 Accuracy of the localized basis as a function of the norm of the 
difference between the global and the localized basis 

We introduce in this subsection lemmas that will be used to control the accuracy of the 
local basis {(j)Y^)i^j\f in approximating solutions of ( |1.1[ ) (i.e., the error associated with 
replacing the global basis {4>i)i^ with a localized version {'^i)i^j\f). Although {ipi)i^j\f 
is arbitrary in the following lemma, it is intended to be a localized version of the global 
basis. 

Lemma 7.4. Let {ipi)i,^ be N elements of'HQ{U). For c G M^, c ^ 0, define 

N N 

f (x) := Cj(/)i(x) and w{x) ■.= ''^^Ciipi{x) (7.8) 

1=1 i=l 

It holds true that 

where the constant C is finite and depends only on Xrainio), Aniax(ffl) CL^d Jl. 



Proof. We have, using (2.31), 



c^Mc 1 

h-M\nl(n) = (^T^)'l|div(aV?;)||^2(n) (7-10) 
where M is the N x N matrix defined by 

Mij := [ iV(l)i - V^ifiVcPj - Vc^j) (7.11) 
Jn 



and P is defined in (2.13). Note that since is symmetric positive definite, there exists 
a symmetric positive definite N x N matrix such that @ = T^. Recalling that P = Q~^, 
we obtain 

c^Mc 

sup = X^UTMT) < Ama.(M)Ainax(e) (7.12) 



Using (2.6) we have for c G 



t,N 



N 

c^Qc= / {^G{x,Xi)ci)^ dx (7.13) 
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It follows from the Cauchy-Schwartz inequality that 

N 



Using (2.11), we conclude that 



Amax(e) < NK 



(7.14) 



(7.15) 



where the constant K is finite and depends only on Amin(a), Amax(a) and 0. To sum- 
marize we have obtained that 



W — ww-iji 



Him 



div(aVt;) 



<CiVl/2(A^ax(M)) 



(7.16) 



IL2{Q) 



We now need to control Amax(-^j'^)- Since M is symmetric and positive its maximum 
eigenvalue is less or equal than its trace. It follows that. 



N 



,{M) Ma < N max Mi 



Combining (7.17) with (7.16) we deduce (7.9). 



(7.17) 



□ 



The following lemma bounds the accuracy of the localized basis {4>^°'^)i£j\f in approx- 



imating solutions of (1.1). In particular, it shows that if maxjg_yv 



iloc I 



IS 



sufficiently small then the accuracy of the localized basis is similar to that of the global 
basis. 



Lemma 7.5. Let u he the solution of equation (1.1), and n^'^"*^ he the finite element 



solution of (1.1) over the linear space spanned hy the elements {(f), 
identified in (7.3). It holds true that 



loc 



.N} 



\U - U^^'^'^W^.^^^ < C\\g, r2 



T-i,r^. 1 H + A^max Wchi 



iloc I 



(7.18) 



where the constant C depends only on Aniin(ffl); Amax(a) CLnd the (uniform) hound on the 
aspect ratios of the elements (tetrahedra for d = 3, triangles for d = 2) ofTn- 



Proof. Let u^^ be the interpolation of the solution of (1.1) over the elements i = 
1, . . . ,N} identified in (2.4). Let be the interpolation of the solution of (1.1) over 

the elements {f/*'"^, i = 1, . • • , N} identified in (7.3). Using the triangle inequality we 
have 

(7.19) 



\n - u'-''°^\\^i^^^ < \\u - 



1 + u" 



„,in,loc II 



Using Corollary 3.2 and Lemma 7.4 we deduce that 



\u-u'-''°'\\Hi^^^<CH\\g\\^2^^^ + CNum^^^^^ (7.20) 
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Using the minimization property given in Theorem |2.6| we have 

II div(aVu-)||L.(^) < II div{aVu)y^^^ = \\g\\^.^^y (7.21) 

We deduce that 



\\n - n-^'''''\\nm < CMj^^^^^ [H + iVmc^x ||^, - 0- j . (7.22) 
We conclude the proof using the fact that u^'^°^ minimizes 

over ah v in the Hnear space spanned by the elements {(/>'°'^, i = 1, ■ ■ • , N}. □ 

7.3 Identification of the difference between the global and the localized 
basis 

In this subsection we will represent 4>i — as the solution of a fourth order PDE. 
We will assume that d^i n O is at some strictly positive distance from the set of coarse 
nodes, i.e. writing 

8i := inf llx — xJI (7.23) 

we assume that > 0. Although the construction of the localized basis 93-°'^ does not 
require 5i to be strictly positive our proof of the accuracy of the localized basis will require 
that 5i is uniformly bounded from below by some arbitrary power of H (i.e., C for 
some C > and k > 1). More precisely, our (final and simplified) a posteriori error 
estimates (Theorem 7.15| ) will be given under the assumption that minjg_/^5j > H^\^/1Q 
where 

Fmin := min - (7.24) 

and that the assumption that maxjg;\/5j ^ H <\ (we also make these assumptions to 
simplify our expressions, note that there is no no loss of generality here since O can be 
re-scaled to have a diameter bounded by 1). 

Let / G H~^{n) and consider the system of equations with unknowns c € and 
X G Vo (where Vq is defined in (2.26)) 

div ^aV(div(aVx))) = / + Zlili ^4(5(2; - Xj) onO. 

X = div(aVx) = on d9. (7.25) 

X = on {xj : j G A^}. 



Proposition 7.6. Equation (7.25) admits a unique solution x G ^ given by 

Xix) = x\x) + x\x) (7.26) 

where 

X\x)= [ T{x,y)f{y)dy (7.27) 
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and 



N N 

X^{x) = - '^'^T{x,Xj)Pj^kX^{Xk). 

k=l j=l 



(7.28) 



Proof. Noting that r is the fundamental solution of (2.15) we obtain that (for c G M^) 
X is the unique solution of 



div \ aV ( div(aVx)) j = / + Yln=i '^i^ (x - Xi) on 0, 
= cliv(aVx) = on dQ 



if and only if 



N 

X{x) = / T{x,y)f{y)dy + Y] 



T[X, Xj)Cj. 



The additional constraint 
is then equivalent to 



X = on {xj}jeAf 



N 



(7.29) 

(7.30) 
(7.31) 
(7.32) 



X^{Xi) + '^T{Xi,Xj)Cj = 

i=i 

which (recalling that T{xi,Xj) = Qij and that P = Q^^) admits the unique solution 
c = —P{x^{xi))i,^. Note that xi has well defined values at the nodes {xj)j^j^ since, by 
[651 133] . it is Holder continuous over Q.. □ 



Consider the equation 

div (^aV(div(aVx)) 
X = div(aVx) = on dO, 
X = on {xj : j G J\f}. 



div ( aV( div(aV0j°^)) ) + J^Zi - ^i) on ^ 



(7.33) 

We will show we will in the proof of Lemma 7.7 that (7.33) admits a unique solution. 
We will denote that solution by Xi arid obtain (in Lemma 7.7) a representation formula 
for it. 



Lemma 7.7. It holds true that 

<t^^ = (t>'r + X^ (7.34) 

where Xi = Xi + Xi with (n is the unit surface vector oriented towards the exterior of 



x]{x) 



dQinn 



and 



r(x,2/)n-oV(div(aV<^J°"))dy- [ G{x,y)n ■ aV (pf^ dy (7.35) 

(7.36) 



N N 

xKx) ■■= - ^^'r{x,Xj)Pj^kx}{xk)- 
k=i j=i 
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Proof. Initial proof. Let / be an arbitrary element of T-l ^{0.) and x the solution of 

(7.37) 



( 7^ and ( |7.3lD . We have (using cpi - (^f^ = on {xj : j £ M}) 

jj,ct>, - <l}r)f = ^(0. - <t^'r) div (aV( div(aVx))) . 



Integrating by parts we deduce that 



0^)/ = / div (aV(</), - cl^n) div(aVx) + / div(aVx) n • aVcj^r. 



n/dUi 



dn,nn 



(7.38) 

Integrating by parts again (using the continuity of div(aV(/)^°'^) across d^i n Q and the 
fact that X is null on the nodes Xj) we have 



[ div (aV((/>i - ^f")) div(aVx) = / x « • aV div (aV(/>l°") . 



(7.39) 



Replacing x by the solution expressed in (7.26), i.e 



N N 

dy -'^'^T{x,Xj)Pj^k j T{xk,y)f{y)dy (7.40) 

k=l j=l 



we obtain that 



G{x,y)f{y)n-aV(t>'r{x)dydx 



yen Jx^dVLiHU 



+ 



yen Jxe9n,nn 



{x, y)f{y) n ■ aV div {aV4>rix)) dy dx 



N N 

k=i j=i 

N N 

fc=i j=i 



-{x, Xj)T{xk,y)f{y) n ■ aV div (^aV(l)°^{x)) dy dx 



(7.41) 



G{x, Xj)T{xk, y)f{y) n • aV^-°'=(x) 



yen JxedUinQ 



which leads us to (7.34). 

Alternate proof. For the sake of clarity we will now also give an intuitive proof 
of ( |7.34 ) based on a direct application of Proposition 



7.6 



Note that since (f>^°'^ and 

div(aV0j°^) are continuous across dQi CiQ we have (using the representation formulae 
( [2T6I ) and that 



div ['aV(div (aV((^i - 0^))) ] = / + ^Q^C 



(7.42) 



i=l 



where 



f = n- aV{diY{aV4>'n)^dn.nn + div (aV(n • aV4>'rSdn,nn)) (7.43) 
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and (5gn,nr2 is the Dirac (surface) distribution on 517 j n 17. We conclude by using 
Proposition |7.6| and integrating by parts. Note that although Proposition |7.6| requires 
/ G ^^^(ri), our initial proof shows that it can still be applied here (with / G ^^'^(17)) 



to show (7.34). Note that our assumption that d^li D O (the support of /) is at some 
strictly positive distance from the coarse nodes {xj)j^^ (i.e., 5i > 0) implies xj has well 
defined values at the coarse nodes {xj)j^j^ (since div(aVXi) is in TH} in a neighborhood 
of those nodes). □ 

7.4 Reverse Poincare inequality 

The following lemma is a generalization of the classical Caccioppoli's inequality (usually 
referred to as the reversed Poincare inequality, see [32] and Proposition 1.4.1 of [51j). 
We will remind its statement and proof for the sake of completeness. 

Lemma 7.8. Let Di,D2 be open subsets of Dq such that D2 C Di (and such that Dq/Di 
is non-empty). If v € 'H}{Dq) satisfies div(aV?;) = in Di, and v = on ODq, then 



^^(^-)- dist(Z).,ZJo/IJi) ''"''^^(^^) + UI^^"(°^")|l^-(^J'"''^^(^^)J • (^-'^^ 

Proof. The idea of the proof is as follows, choose r/ to be a differentiable function on Dq 
such that ?7 = 1 on 1)2, 7? = on Dq/Di, < r/ < 1 and ||Vr?||ioo(OQ) < distCDa^o/J'i) ' 
Using the fact that rfv = on dDi we obtain that 

/ V{jfv)aVv = -i rfv<X\\{aVv) (7.45) 



which leads to 



/ rf'VvaVv = —2 I rjvVrjaVv — / ?7^v div(aVv). (7.46) 
Jdi Jdi Jdi 

Using the Cauchy-Schwartz inequality, the uniform bound on Vr/, and the uniform bound 
on a we deduce that 

/ Tj'^VvaVv < , / n '^n /n ^ ( / ^^^^«^^) ' II^IIl2(z)i) 
Jdi dist(L>2,L>o/^i) Jdi ^ ' (7.47) 

+ II div(aV7;)||^2(^^^||t;||i2(£,^). 

1 

Writing ( J^^ rj^VvaVv^ ^ , we conclude by using the fact that, for positive x, y,z, x < y+ 
z/x implies that x < y+\/z and also using the fact that J^^ VvaVv < rfVvaVv. □ 

7.5 Bound on the maximum eigenvalue of P 

In this subsection we will provide an upper bound on the maximum eigenvalue of P as 



defined in (2.13). 
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Proposition 7.9. Let Q be defined as in (2.6), P as in (2.13), and //mm o-s in (7.24). 
It holds true that 

N 

(7.48) 



Amax(-P) < 



Proof. The definition (2.13) implies that the maximum eigenvalue of P is the inverse of 
the minimum eie envalue of 6 as defined in Ko^. Let / G such that \\l\\ = 1. We 
wish to bound from below l^Ql. Write 



N 



){x) := '^G{x,Xi)lj 



(7.49) 
(7.50) 



and observe that 

fei = \\v\\h^ay 

Since ||/|| = 1, there exists i ^ Af such that > Let Bp be the intersection of Q 

with the ball of center Xi and radius p. Let p = i/min/8. Observing that v is harmonic 
in Bjp/Bp we have using Caccioppoli's inequality (see lemma 7.8) that 



\^v\\L2(B5p/Bip) < 



Using Green's identity and (7.49) we have 



n ■ aVv 



OBbp 



PIIL2(B6p/B3p)- 



div(aV?;) 



(7.51) 



(7.52) 



Let ?7 be a smooth function equal to one on dB^p, zero on B^p and such that < < 1 
and II Vr/||x,oo(Bg^) < C/p. Integrating by parts, we have 



n ■ aVv 



dB. 



5p 



dB^p 



B^pl B4p 



(7.53) 



VrjaSIv. 



Observing that 



B^p/Bap 



Vr]aVv\ < - II 1^2(^5^/5^^) 



P 



we obtain from |/j| > l/\/N, (7.52) and (7.53) that 

\l2{Bsp/B4p)- 



1 C',,^ , 

< — 



N p 



Combining (7.55) with (7.51) we deduce that 

^2 



C 



P 



N 



< \\v\ 



LHBep/Bsp)- 



(7.54) 



(7.55) 



(7.56) 
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Combining (7.56) with (7.50) we deduce that 

Fei > (7.57) 
which concludes the proof. □ 

7.6 Pointwise estimates on solutions of elliptic equations with discon- 
tinuous coefficients 

Write B{x, p) the Euclidean ball of center x and radius p and Q{x, p) the intersection of 
B{x, p) with ri. 

Lemma 7.10. Let v € Ti^lQ,), d < 3 and xq G Then there exists a constant C 
depending only on Ainin(fl) o-n-d Aniax(o) such that 

C f ^ d 

H JQ{xQ,2p) 



Proof. Lemma 7.10 is classical and we refer to [66] for its proof. The idea is to decompose 
V as f = vi + f 2 where vi and V2 satisfy div(aVwi) = div(aVu) and div(aVt>2) = on 
Q{xo, p) with boundary conditions t;i = and V2 = v on dQ{xo,p). By the estimates on 



the L^-norm of the Green's functions given in Lemma 2.2 (see also |66| 1651 [33]) 

II^^i||l-(o(xo,p)) < Cp^-^\ div{aVv)\\^,^^^^^^^^^y (7.59) 

Since V2 is harmonic we have (using for instance Theorem 5.1 of [SS]), 

\\v2\\L^{n{xo,p)) ^ C'lb2||L2(n(xo,2p))- (7.60) 

□ 

7.7 Control of the norm of the difference between the global and the 
localized basis 

Write dist(x,^) the distance from a point a; to a set A in Euclidean norm. Define for 

nf ■={xeni\ dist(x, dQi nQ) <6}. (7.61) 
Lemma 7.11. We have for 5i < 1 

\\xl\\Lo.(n)<C6i^dw{aV<p'n\\ +C6;^\<p'r\\ 3.,. (7.62) 

Proof. Let r/ be a smooth function equal to 1 on Q-^ , on ilj / and such that < rj < 1 
and such that 

C 

sup |Vr/(x)| < - (7.63) 
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for some constant C independent from H and Jj. Let be the function defined in 



(7.35). Noting that r/ = 1 on d^i H Q, T{x,y) = and G{x,y) = for y G 90 we have 
r(x,y)n-r?aV(div(aV(/>l°"))dy- / G(x, y)n • r/aV(/>l°" dy. (7.64) 

J dQ.i 



We deduce, after integration by parts and using the fact that div ( aV( div(aV(/); 



on Q-^ , that 



■■ I Vr(x,y)??aV(div(aV(/)l°^)) dy+ f r(x, y)V7?aV( div(aV0j°^)) 



Hi 



iloc 



(7.65) 



G(x,y)r/div(aV0^°=)dy. 



Noting (by integration by parts) that 



VG(x, y)r?aV0i°^ dy = -<pr{x)r^{x) + / VG(x, y)aV?7<Ar ^^2/ 



/ loc 



iloc 



<|</.l°^(x)r/(x)| + C||Vr/,/.M| ||VG(x,y)|| 



(7.66) 



and observing that the norms of VT(x,y), T{x,y), G{x,y) and the norm of 
VG(x,y) are bounded by a constant G depending only on Amin(a)5 Amax(a) and 17 (see 
proof of Lemma 2.2, equation ( |2.14 ), |39| I67j and references therein) we deduce that 

V(div(aV<Al°^)) 

G||div(aV(/)J°^^" 



+ G5ri||V(/.l°^' 



L2{Q^2 ) 



III ^ 

3c5,- 



'1 II j,loc| 



(7.67) 



Since (div aV0j°'^) is harmonic on il-" we obtain from Caccioppoh's inequahty (i.e. 
with Do = J^i, L»2 = and Di = ft-^ ) that 



Lemma 



7.8 



V(div(aV<Ar)) 



<G(^ri||div(aV(/>l°=)|| 3., . 



(7.68) 



Therefore 



\X' 



.ll^oc.n)<C<5r'l|div(aV0i°^)|| a., 

+ G5ri||v<Al°^|| +G<5, 



-lIUloci 



(7.69) 



L°o(a2 ) 
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Similarly, using the reverse Poincare inequality (Lemma 



35,; 



7.8 



with Dq = Qi, D 



i, ^2 



and Di = ) and Young's inequality (|o6| < + -6^) we obtain that 



loc I 



< C6: 



Therefore 



X: 



,111 

* llL°=(f7) 



-111 j.loc| 



<C(5r2||div(aV0; 



35,- 



+ SAl div(aV0; 



loc\ 



35^ . 



loc\ 



+ C6, 



-211 -loci 



3c5j 



-HI j,loc| 



(7.70) 



(7.71) 



Now to control 



^ote that on Q,-^ , (f)f"^ can be decomposed as 

iloc\ 



iloc 



vi + V2 where vi and V2 satisfy div(aVwi) = 6i\{aV<j)^) and div(aVf2) = on Vt^ with 

1 ^ 

6:-°'^ on dVt.^ . We therefore obtain that 



boundary conditions vi = and V2 

1 1 j,loc II 



L°°{nf) L°^{nJ) L°°(n7") 



(7.72) 



By | .65t [33] (one can also use the bounds on the Green's function given in the proof 
Lemma |2.2|) we have 



bill h < C||div(aV(/>l°^)|| . 



Since V2 is harmonic we have, 



iloc I 



Using lemma 



7.10 



iloc I 



Il0i llL°°(B(a;o,p)nQi) - V 



1^211 <c\wi-\\ 



we obtain that for xq G H i7j and p = 6i/8, 
C 



(7.73) 



(7.74) 



< 



P JB{xo,2p)r\ni 



^-i"')^) ^ + II div(aV(/>l°=)||^2 



{B{xo,p)nni)- 
(7.75) 



Combining (7.75) with (7.74) we deduce that 



1^^211 <C6. 



iIoc I 



, +S; ^||div(aV(/>^) 



Using (7.72) we deduce that (for 6i < 1 and d < 3) 

d 



i loc I 



s <C5, ^ 



jloc I 



+ C div(aV(/>; 



loc\ 



3Si . 



Combining (7.77) and (7.71) we deduce (7.62) for 6i < 1. 



(7.76) 

(7.77) 

□ 
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Lemma 7.12. Let xf be defined as in (7.36). We have 



W'^xlWmn) < N\\xi\ 



max \/ Pjj . 



Proof. We have 
Using 

we deduce that 



N 



k=l 

W^Mmn) < C||div(aV0fc)||^2(Q) 

N 



k=l 

Lemma 7.13. We have for < (5j < 1 

W^xhlmn) < C6r^{\\div{ay^n\ 

Proof. Using 



iloc I 



VvaVv = — t>div(aVf 



we obtain from (7.35) that 



with 



and 



VxJaVx! <h + l2 



h= I / xi{x)G{x,y)n-aV{dW{aV^'n)dy 



xi{y)n-aV4'^r dy. 



dn.nn 



Defining rj as in the proof of Lemma 7.11| we have (since rjxi = on 90 j) 



Integrating by parts, we obtain 



Xi(y)r/n • aVcpf"" dy 



h 



rixi div(aV0j°'=) - / riVxiaV^\ 



loc 



XiVvaVct^r- 



Using Young's inequahty we obtain 



/ r,VxiaVcl}r\<\ [ VxiaVxi + C\\V4>'rf 



(7.78) 

(7.79) 
(7.80) 

(7.81) 

□ 

(7.82) 

(7.83) 

(7.84) 
(7.85) 
(7.86) 

(7.87) 
(7.88) 

(7.89) 
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Therefore 



I/2I <||xi||L-{Q)||div(aV0i°^)|| +C5-i||xi||L->(n)l|V(/>M| 



+ C\\V<l}rf +] [ VxiaVxi 

L2(n7) ^ Jn 



(7.90) 



Similarly we have 

I , 

Jxen Jyedfi 

which, after integration by parts, leads us to 



Xiix)G{x,y)r]n- aV( div(aV</>-°'=)) dy 



(7.91) 



Xi[x 



)G{x,y)vdw (aV{div{aV4>'n)) dy 



+ 



+ 



I x£f2 J ydili 



Xi(x)G(x,y)Vr?oV(div(aV(/>l°=)) dy 
Xi{x)VyG{x, y)7]aSl{diY{aV(t>^r)) dy. 



(7.92) 



Using the fact that rj div ( aV( div(aV0'°'^)) j = (since the support of t/ does not contain 



any coarse nodes), and using the bounds the L^-norm of G (see proof of Lemma 2.2) as 
well as 



( / Xi{x)yyG{x,y)f < C||xi||ioo(f2) 



we obtain that 



h<C{l + 6r^)\\xi\\L^in) aV( div(aV</.J°^)) 



(7.93) 



(7.94) 



Combining (7.84), (7.90) and (7.94), we have obtained that 

[ VxjaVx} <G\\xi\\L-(n)\\div{aV4>'n\\ +G\\Vct>'rf 



+ C(l + <5ri)||xi|U^(f,)||v(div(aV<Al°^)) 
+ C6r^\\xi\\L^in)\\V<i^°'" 



(7.95) 



L2(nf ) 

Therefore, using Young's inequality, we have obtained that 
pxlWmu) <C^(l + 5r')llxi||L-(n) + ||div(aV(Al°^)|| +C||V0 

V(div(aV(/>l°^)) 



loc I 



+ 



We deduce ( |7.82D by combining ( |7.62D and ( |7.68D with ((7J0). 



(7.96) 



□ 
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Combining lemmas 7.7, 7.12 and 7.13 with Proposition 7.9 we deduce the following 
lemma. 



Lemma 7.14. For < Si < 1, it holds true that 
||V(<A^ - '^1°^)|L.(^) < CN'H^J-HW div(aV,^^)|| ss, + \\cp'r\\ ). (7.97) 



7.8 A posteriori error estimates. 



Theorem 7.15. Let u be the solution of equation (1.1), and u ' be the finite element 



solution of (1.1) over the linear space spanned by the elements {<i)i, i = 1,...,A^} 



identified in (7.3). Assume that vnvni^j^ 6i > -fTmin/lO (where 6i is defined in (7.23) and 



Hmin in (7.24)^ and that maxjg_/v5j < H < 1. It holds true that 



u 



HAoc I 



{7: 



where the constant C depends only on \rmn{o), Amax(a) (^nd the (uniform) bound on the 
aspect ratios of the elements (tetrahedra for d = 3, triangles for d = 2) ofTn and 



E = H. 



-7-3d 



maxi^\\div{aV4>Yl\\L2^nH) + 



I loc I 



(7.99) 



where Q,^ is defined in (7.61). 



Proof. Theorem |7.15| is a straightforward consequence of Lemma |7.14| and Lemma |7.5 

mm 



and the inequality < CH^- . 



□ 



Remark 7.16. Observe that there is no loss of generality in requiring that mmii^j\f 6i > 
-f^min/10 and that maxjg_/^5j < -ff < 1- Note also that minjg_^5j > f^min/lO corresponds 
to the condition that the minimum distance between dQi D 17 and the set of coarse nodes 
is at least i/min/lO. Although we need this assumption for the simplicity of the proof, 
numerical experiments suggest that it is not needed for convergence with an optimal rate. 



Remark 7.17. Numerical experiments show that div(aVi;^-°'^) 



+ 



decays exponentially fast as a function of Ui where Ui is the minimum number of (coarse) 
edges (i.e., number of layers of triangles) separating Xi from dVLi n Hence if Ui is 
of the order ofln{l/H) we can derive a posteriori error estimates from the posterior 
observation that E < H . We will in a sequel work derive a priori estimates by bounding 
the exponential decay o/|| div(aVi?!)^°'^)||^2(f^H-) + ||</''°'^||£^2(f^H) a-s a function of the number 
of layers of triangles surrounding each node Xj. These a priori error estimates will 
guarantee \\u — ^^'^°'^||-Hi(r2) — ^WdWlJ^i^^)^ provided that minjg_/^ dist(a;j, (9Qj n O) > 
C*H\x\^ for some constant C* depending on Ainin(o) o,nd Amax(ffl)- We refer to this 
property as super- localization. Note that each (super-localized) basis function 0^°'^ is the 
solution of a sparse/banded/nearly diagonal linear system. This property is essential in 



keeping the computational cost minimal (see Subsection 10.3). 
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Figure 1: Mesh details. 




(a) Primal and dual (b) Xi and Qi with (c) Qi with I = 2 (d) Qi with I = 3 
mesh I = 1 



8 On time dependent problems 



As shown in sections 4 and 5 of |56| the accuracy of the basis elements and 
remains unchanged when those elements are used to approximate the solutions of the 
parabolic and hyperbolic equations associated with — div(aV), i.e., for equations of the 
form 



and 



jdtu{x,t) — div (^a{x)'\7u{x,t)j = g{x,t) {x,t) G Qt] 9 G L'^i^x), 

= on d^^T: 

p{x)dtu{x,t) - div (^a{x)Vu{x,t)^ = g{x,t) {x,t) £ Qt] 9 G L'^i^x), 

u = on O^It, 

dtu = on nx{t = 0}. 



More precisely the solutions of (8.1) and (8.2) are approximated, via finite element 
methods, using the linear space spanned by the interpolation basis as test and 
solution space. Note that the approximate solutions are of the form 

u^(x,t) = (8.3) 



We refer to sections 4 and 5 of [56] for precise statements and proofs. 



9 Numerical implementation and experiments 
9.1 Implementation 

Our method is, in its formal description, meshless in the sense that only the locations 
of the points {xj)j^j\f are needed to construct the rough polyharmonic splines 0j or 
as solutions of the minimization problem (2.4) or (7.3). 
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For numerical implementation, we need to formulate and solve the problem in the 
full discrete setting. Although, for the sake of conciseness we have provided an error 
analysis only in the continuum setting, this analysis can be extended to the fully discrete 
setting without major difficulties. 

We will now describe the implementation of our method in the discrete setting. Let 
Th be a tessellation of of resolution H with nodes {xj}. Assume that h H and that 
Th is obtained by repeated refinements of Th- 

Let Wh{^) be the set of piecewise linear functions on Th with zero Dirichlet boundary 
condition on dQ. For E W/i(0), we can define gu, which approximates — divaVw on 
Jl. gu is piecewise constant on the dual mesh of Th, and can be defined by the following 
finite volume formulation, 

guiVi) = j^[ Vlv.aVu (9.1) 



where Vj is the dual Voronoi cell associated with the node of Th and ly, is the 
characteristic function which equals to one on Vi and zero elsewhere. 



Subfigure 2(a) shows the primal (finite element) mesh and its dual mesh used for the 



numerical implementation in this paper: the red dots are the nodes Xi of the triangulation 



Th- In Subfigure 2(a), the red square around each yi is the Voronoi cell Vi associated 
with the node yi, and the dual mesh is composed of such cells. 
With this definition, the local basis 



tional formulation (9.2) 



can be readily computed through the varia- 



Minimize 
Subject to 




4> G Wh{n,),(t> = on dni u n\ni 

= 5ij,jGM. 



(9.2) 



9.2 A two dimensional example. 

In this example we consider = (0, 1) x (0, 1) and 

1 l.l + sin(27rx/ei) 1.1 + sin(27ry/e2) 1.1 + cos(27rx/e3) 

6 1.1 + sin(27ry/ei) 1.1 + cos(27r2;/e2) 1.1 + sin(27ry/e3) 

l.l + sin(2W.4) ^ l.l+cos(2W..) ^ ^ 
1.1 + cos(27rx/e4) 1.1 + sm(27ry/e5) 

where ei = \,e2 = ^,£3 = 3^,64 =^,£5 = l^- 



(9.3) 



Note that a, as defined by (9.3), corresponds to the one used in [55] (Example 1 
of Section 3) and |50j. The number of coarse nodes is N = Nc x Nc with Nc = 32. 
The (regular) coarse mesh Th = T^ is obtained by first (uniformly) subdividing Q, into 
Nc X Nc rectangles, then partitioning each rectangle into two triangles along the direction 
(1,1). We can further refine the mesh to obtain T^, T^, T^. Global equations are 
solved on a fine triangulation with 66049 nodes and 131072 triangles. We compute 
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localized elements on localized sub-domains il' defined by adding I layers of coarse 
triangles around Xj. More precisely ilj is the union of triangles sharing Xj as a node and 
is the union of Q[ with coarse triangles sharing a node with a triangle contained in 
O'. We refer to subfigures 2(b) , 2(c) and 2(d) for an illustration of with I = 1,2, 3. 



We refer to Figure [2] for the results of our numerical experiment. Subfigure 3(a)| 
shows (f)i (one element of the global basis in dimension two). Subfigure |3(b) shows 
a slice of (pi along the x-axis. Subfigure 3(c) 



shows 



aIoc I 



m 
1,. 



ri^, and L° 



norms in the log scale as a function of the number of layers I = 1, ... ,8 used in the 
localization. Note that ||(/>j — (/'^°'^|| decays exponentially fast with / and that the decay 
rate is nearly independent from the choice of the norm. Subfigurea |3(d)| and |3(e) show 



\u — u 



H,loc I 



\l2(q) and II M — u^'^°'^\\H^(n) in log-scale as a function of the number of layers 
/ = 1 , . . . , 8 (defining the size of localized sub-domains) for different values of (coarse) 
mesh resolution H = 1/4, 1/8, 1/16, 1/32. Note that as the number of layers increases. 



u — u 



H,loc\ 



decreases first and then quickly saturates around 0{H). 



9.3 Wave equation 



In this example we consider a and Q as defined in Subsection 9.2 and we compute the 



solutions of wave equation (8.2) up to time T = 1. The initial condition is u{x,0) = 
and ut{x,0) = 0. The boundary condition is u{x,t) = 0, for x £ dQ. The density p 
is uniformly equal to one and we choose g = sin(7rx) sin(7ry). Subfigures [4(^ and [4(b)] 



show u{x, T) and u^'^°^{x, T) with T = 1 (n^''°'^(x, T) is computed over the (approximate 
solution) space spanned by the localized basis elements (/^"^ obtained by choosing Qi as 
the union of 3 layers of triangles around Xi (hence, we choose / = 3). The resolution of 
the fine mesh is h = 1/80. The resolution of the coarse mesh is H = 1/10. Subfigure 



4(c) shows approximation error \\u — u ' (.,r)|| in L , Ti , and L°° norms in the log 



scale as a function of the number of layers 1 = 1, 



9.4 A one dimensional example 

In this example we consider = (0,1) and 

1 ^ 

a{x) := 1 + - sin ( ^ fc "(Cifc sm{kx) + C2k cos(fcx))) (9.4) 
fe=i 

where {Cifc} and {C2fc} are two independent random vectors with independent entries 
uniformly distributed in [— 5, 5]- This example is taken from \A'2\ I5U| and it has no scale 
separation. Note that a has an algebraically decaying (Fourier energy) spectrum, i.e., 

{\a{k)\')^\kr, (9.5) 

where ^-^ denotes ensemble averaging, and a{k) is the Fourier transform of the coefficients 
a(x). We choose = 20 in the numerical experiment. 
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Figure 3: Wave equation of Subsection 



9.3 
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Figure 4: One dimensional example (9.4). 



local basis at node 40 



log |10 +|d|) alnode 40 



0.1 0.2 0.3 OA 0.5 0.6 0.7 0.8 0.9 




0.1 0.2 0.3 OA 0.5 0.6 0.7 0.8 0.9 



(a) 



(b) (jii in log scale 



)r of local basis at node 40 



0.4 0.6 



a Qiad of basis al node 40 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 



(c) (|)^ 



(d) -div(aV(/)i) 



logi of -div a grad ol basis al node 40 





(e) — div(aV<^i) in log scale 



(f ) Matrix Q 




Matrix [og(1+abs(P)) 



(g) Matrix P 
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(h) Matrix P in log scale 



The coarse nodes {xj)j^j\f correspond to = dc = 80 points uniformly distributed in 
(0, 1). More-precisely, those points correspond to the nodes of a coarse mesh of (0, 1) and 
the distance between two successive points is H = 1/81. The fine mesh is then obtained 
by refining (bisecting) the coarse mesh 3 times. Local bases are computed on intervals 
with I layers, where I = 1, . . . , 8. More precisely our localized sub-domains (over which 
the elements (/)'°^ are computed) are Qi = (xj — //81, Xj + //81) n for / = 1, . . . , 8. 

We refer to Figure |4] for the results of our numerical implementation. Subfigure 



5(a) shows the local basis (1) °^ for the node i = 40 for various degrees of localization 
(i.e., for I = 1,...,8). Subfigure |5(b) shows the global basis in the log scale (i.e., 
logjn(10~^+|(/)4o|)) and illustrates the exponential decay of </>i(x) away from Xj. Subfigure 



5(c) 



shows (f>i — (j)^'^ , the difference between the global and the local basis, at node i = 40, 
for various degrees of localization (i.e., for / = 1, . . . , 8), and illustrates the decay of this 
difference with respect to the number of layers / defining Q 
— div(aV(/)j) for i 



40. Subfigure 5(e) shows 



(i.e., log;^o(10 +1 divaV0j|)). Subfigures 5(d) and 5(e) illustrate the exponential decay 



of — div(aV(/>j)(x) away from Xj. Subfigure 5(f) shows the matrix Qij, as defined in 



Subfigure 5(d) shows 
div(aV0i) for z = 40 in the log scale 



(2.6), as a function of (two dimensional surface plot). Similarly Subfigure 5(g)| 



shows the matrix Pij, as defined in (2.13), as a function of Subfigure 5(h) shows 



P in the log scale and illustrates the exponential decay of the entries of P away from 
the diagonal. 



Analogy with polyharmonic splines. Note that as H the matrix Q shown in 



Subfigure |5(f) converges towards T{x,y) the fundamental solution of the operator L = 
(div(aV-)) . Note also that since P^^ = Q, P is, in fact (in analogy with [591 [601 EI]), 
a discrete approximation of the operator L (i.e. X^jLi Pi,ju{xj) can thought of as an 
approxim ation of Lu{xi)). Note the localization of P (near its diagonal) as shown in 



Subfigure 
the discre 
elements c 



5(g) 



Note also since 



N 



1 ^ij^i.'^: ■^j 



is obtained by applying 



;e version of L against the fundamental solution of L and in that sense the 



or 



tloc 



see Subfigure 5(a) ) are approximations of masses of Diracs. 



10 On numerical homogenization 
10.1 Short overview 

As mentioned in the introduction the field of numerical homogenization concerns the nu- 



merical approximation of the solution space of ( 1.1 ) with a finite-dimensional space. This 



problem is motivated by the fact that standard methods (such as finite-element method 
with piecewise linear elements [10]) can perform arbitrarily badly for PDEs with rough 



coefficients such as (1.1). Although some numerical homogenization methods (described 



below) are directly inspired from classical homogenization concepts such as periodic ho- 
mogenization and scale separation [12], oscillating test functions, G or i/-convergence 
and compensated compactness [52l EH |33] and localized cell problems fST] , one of the 
main objectives of numerical homogenization is to achieve a numerical approximation of 



35 



the solution space of (1.1) with arbitrary rough coefficients (i.e., in particular, without 
the assumptions found in classical homogenization, such as scale separation, ergodicity 
at fine scales and e-sequences of operators). Numerical homogenization methods differ 
in (computational) cost and assumptions required for accuracy and by now, the field 
of numerical homogenization has become large enough that it is not possible to give a 
complete review in this short paper. For the sake of completeness, we will quote below 
the the non exhaustive list of numerical homogenization methods discussed in [56J : 

- The multi-scale finite element method \4:2\ \70\ HH |27] can be seen as a numeri- 
cal generalization of this idea of oscillating test functions found in ff-convergence. A 
convergence analysis for periodic media revealed a resonance error introduced by the 
microscopic boundary condition [42^ Wl\ . An over-sampling technique was proposed to 
reduce the resonance error |42j . 

- Harmonic coordinates play an important role in various homogenization approaches, 
both theoretical and numerical. These coordinates were introduced in [44j in the con- 
text of random homogenization. Next, harmonic coordinates have been used in one- 
dimensional and quasi-one-dimensional divergence form elliptic problems [HI [7] , allowing 
for efficient finite dimensional approximations. The connection of these coordinates with 
classical homogenization is made explicit in [3] in the context of multi-scale finite ele- 
ment methods. The idea of using particular solutions in numerical homogenization to 
approximate the solution space of (1.1) appears to have been first proposed in reservoir 
modeling in the 1980s [16] , |69j (in which a global scale- up method was introduced based 
on generic fiow solutions (i.e., fiows calculated from generic boundary conditions)). Its 
rigorous mathematical analysis has been done only recently [55] and is based on the fact 
that solutions are ?^^-regular with respect to harmonic coordinates (recall that they are 
^^-regular with respect to Euclidean coordinates). The main message here is that if 
the right hand side of (1.1) is in L^, then solutions can be approximated at small scales 
(in "H^-norm) by linear combinations of d (linearly independent) particular solutions {d 
being the dimension of the space). In that sense, harmonic coordinates are (simply) 
good candidates for being d linearly independent particular solutions. The idea of a 
global change of coordinates analogous to harmonic coordinates has been implemented 
numerically in order to up-scale porous media fiows [261 EH HE] . We refer, in particu- 
lar, to a recent review article [IB] for an overview of some main challenges in reservoir 
modeling and a description of global scale-up strategies based on generic fiows. 

- In |l23l [29] , the structure of the medium is numerically decomposed into a micro- 
scale and a macro-scale (meso-scale) and solutions of cell problems are computed on the 
micro-scale, providing local homogenized matrices that are transferred (up-scaled) to 
the macro-scale grid. See also [3 [301 1]- We also refer to |3S! for convergence results on 
the Heterogeneous Multiscale Method in the framework of G and F-convergence. 

- More recent work includes an adaptive projection based method which is 
consistent with homogenization when there is scale separation, leading to adaptive al- 
gorithms for solving problems with no clear scale separation; finite difference approxi- 
mations of fully nonlinear, uniformly elliptic PDEs with Lipschitz continuous viscosity 
solutions [18] and operator splitting methods [Sj H] . 
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- We refer to jT5l [H] (and references therein) for recent results on homogenization of 
scalar divergence-form elliptic operators with stochastic coefficients. Here, the stochastic 
coefficients a{x/£,uj) are obtained from stochastic deformations (using random diffeo- 
morphisms) of the periodic and stationary ergodic setting. 

- We refer to [11] for a detailed study of the fluctuation errors of MsFEM and HMM 
in the context of stochastic homogenization and an a rigorous assessment the behavior 
of those multiscale algorithms in non-periodic situations. 



10.2 On Localization. 

The numerical complexity of numerical homogenization is related to the size of the sup- 



port of the basis elements used to approximate the solution space of (1.1). The possibility 
to compute such bases on localized sub-domains of the global domain without loss of 
accuracy is therefore a problem of practical importance. We refer to [19], |24] . [8], [56] 
and |l8] for recent localization results for divergence-form elliptic PDEs. The strategy 
of [19j is to construct triangulations and finite element bases that are adapted to the 
shape of high conductivity inclusions via coefficient dependent boundary conditions for 
the subgrid problems (assuming a to be piecewise constant and the number of inclusions 
to be bounded). The strategy of [24j is to solve local eigenvalue problems, observing 
that only a few eigenvectors are sufficient to obtain a good pre-conditioner. Both [19] 
and [24J require specific assumptions on the morphology and number of inclusions. The 
idea of the strategy is to observe that if a is piecewise constant and the number of in- 
clusions is bounded, then u is locally Ti^ away from the interfaces of the inclusions. The 
inclusions can then be taken care of by adapting the mesh and the boundary values of 
localized problems or by observing that those inclusions will affect only a finite number 
of eigenvectors. 

The strategy of [8] is to construct Generalized Finite Elements by partitioning the 
computational domain into to a collection of preselected subsets and compute optimal 
local bases (using the concept of n- widths [58j) for the approximation of harmonic func- 
tions. Local bases are constructed by solving local eigenvalue problems (corresponding 
to computing eigenvectors of P*P, where P is the restriction of a-harmonic functions 
from uj* onto u G co* , P* is the adjoint of P, and cj is a sub-domain of surrounded 
by a larger sub-domain oj*). The method proposed in [S| achieves a near exponential 
convergence rate (in the number of pre-computed bases functions) for harmonic func- 
tions. Non-zero right hand sides (denoted by g) are then taken care of by finding (for 
each different g) particular solutions on preselected subsets with a constant Neumann 
boundary condition (determined according to the consistency condition). The near ex- 
ponential rate of convergence observed in [8] is explained by the fact that the source 
space considered in [8] is more regular than (see [8]) and requires the computation 
of particular (local) solutions for each right hand side g and each non-zero boundary 
condition, the basis obtained in |8] is in fact adapted to a-harmonic functions away from 
the boundary). 

The strategy of [56] is to use the transfer property of the flux-norm introduced in 
[T^ (and the strong compactness of the solution space |Sn]) to identify the global basis 
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and then (inspired by an idea stemming from classical homogenization [361 1371 1571 17T] ) 
localize the computation of the basis to sub-domains of size 0{^/H)ln{l/H) (without 
loss of accuracy) by replacing the operator — div(aV-) with the operator ^ • — div(oV-) 
in the computation of this basis with a well chosen T balancing the created exponential 
decay in the Green's function with the deterioration of the transfer property. 

In [35], it is shown that for scalar elliptic problems with coefficients, there exists a 
local (generalized) finite element basis (AL basis) with 0((log ^)'^^^) basis functions per 
nodal point, achieving the 0{H) convergence rate. Although the construction provided 
in |38j involves solving global problems with specific right hand sides, its theoretical 
result support the possibility of constructing localized bases. 

The strategy of [l8] is to introduce a modified Clement interpolation Th (in which 
the interpolation of u is not based on its nodal values but on on volume averages around 
nodes), and writing the kernel ofln (i.e., the set of functions u such that Ihu = 0), 
identify Vff^ as the orthogonal complement of with respect to the scalar product de- 
fined by a(u, v) = VuaVv. The finite element basis (pi is then identified by projecting 
nodal piecewise linear elements onto Ih- The work [38] shows that the computation of 
can be localized to sub-domains of size 0{H In jj) without loss of accuracy. The price to 
pay in the method proposed in [48j lies in the fact that the identification of the space 
and the related projections lead to dense linear systems (even after localization) whereas 
the method proposed in this paper leads to sparse (banded) linear systems (even after 
localization). We refer to subsection [10.3 for a comparison of the computational costs 
of the methods proposed here, in [56j and [48]. 



10.3 Computational Cost 

Write dif the number nodes of the fine mesh contained in the local domain i7j [dij is the 
number of nodes/degrees of freedom of 7^ H ilj). Write df the number of interior nodes 
of the fine mesh Th and dc the number of interior nodes of the coarse mesh Th ■ Then the 
cost for solving local basis introduced in this paper, is 0{difdc) (because each (f>^°'^ 
is the solution of a banded/nearly diagonal linear system), and the cost for constructing 
the global stiffness matrix is 0{dif{difdc/df)dc) = 0{dfjd^/ df). Notice that dijdc/dj 
is the number of local domains which overlap with fi,. Therefore, the overall cost is 
0{difdc + dfjd'^/df) = 0{difdc{l + difdc/df)), and the dominant part is O^df^d^/df). 

The method in |48j needs to resolve the null space of the modified Clement inter- 
polation operator on local domains; therefore, solving the local basis requires 0{dffdc) 
operations (because local basis functions are solutions of dense linear systems). The 
overall cost will be 0{dfjdc + dfjd'^/dj) = 0{dfjdc{l + dc/df)), and the dominant part 

is 0{dffdc). 

For quasi- uniform triangulations Th and Th, we have ~ and dc ~ H^'^. For 
the method introduced in this paper, to achieve 0(H) accuracy in norm, we need 
dif (Hlog{l/H))~^. Therefore the overall cost scales Ci ~ ^ i°g(V-^) Similarly, the 
localization method in [56] requires dif ~ {VHlog{l/H))-'^ for 0{H) accuracy (because 
each local basis function is the solution of a banded/nearly diagonal linear system), and 
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the overall cost is C2 ~ ( hh contrast, for the method in [H], the overall cost 

is c, ~ C'^^Hf^Y- 

Therefore, we have the following conclusions, 

* 5i ~ ^ Hiog^i/H) )'^' therefore, the method in this paper is superior to the method 
in in terms of cost (for the same accuracy). 

• When h ~ if^, we have C2 ~ C3. Namely, the method in [56], though with a 
larger localization domain, has the same order of cost as the method in |48| (for 
the same accuracy). 
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